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We calculate the optical spectra of silicon and germanium in the adiabatic time-dependent den- 
sity functional formalism, making use of kinetic energy density-dependent (meta-GGA) exchange- 
correlation functionals. We find excellent agreement between theory and experiment. The success 
of the theory on this notoriously difficult problem is traced to the fact that the exchange-correlation 
kernel of meta-GGA supports a singularity of the form a/q 2 (where q is the wave- vector and a is a 
constant), whereas previously employed approximations (e.g. local density and generalized gradient 
approximations) do not. Thus, the use of the adiabatic meta-GGA opens a new path for handling 
the extreme non-locality of the time-dependent exchange-correlation potential in solid-state systems. 



The first-principle calculation of the optical properties 
of semiconductors is a classic and practically important 
problem in electronic structure theory. The difficulty 
stems largely from the critical role played by electron- 
electron interactions, particularly the so-called excitonic 
effects, i.e. the interaction of an electron in the conduc- 
tion band with the hole left behind in the valence band. 
Early calculations [l[ based on diagrammatic many-body 
theory achieved good agreement with the experiment at 
the price of much computational effort. In recent years, 
the problem has been tackled by several authors 0] , who 
made use of state-of-the-art methods such as the GW ap- 
proximation for the electron self-energy and the Bcthc- 
Salpetcr equation for the electron-hole interaction. These 
methods are computationally demanding and not so eas- 
ily adaptable to an emerging new generation of electronic 
materials, e.g. organic semiconductors and long polymer 
chains. 

A promissing alternative to the traditional many-body 
approach is provided by the time-dependent density func- 
tional theory (TDDFT) Q. This approach directly 
targets the density-density (or, in some versions, the 
current-current response function of a fictitious non- 
interacting system, the so-called Kohn-Sham (KS) sys- 
tem, which is so designed as to produce (at least in princi- 
ple) the same density /current response as the physical in- 
teracting system. The elimination of interactions greatly 
reduces the computational effort, but the complexity of 
the many-body problem eventually resurfaces, since the 
quality of the results is crucially determined by the qual- 
ity of the (approximate) exchange-correlation (xc) poten- 
tial v xc (r,t) in which the fictitious non-interacting elec- 
trons move. 

Successes and failures of the TDDFT approach to the 
calculation of optical spectra of semiconductor are well- 
documented. The first difficulty, which has been known 
since the early 1980s, is that the basic local-density ap- 
proximation (LDA) and its semilocal extensions severely 



underestimate the band gap. The problem with the 
KS band gap can be corrected by the use of orbital- 
dependent functionals or the TDDFT approach 
can be implemented on top of a band-structure obtained 
by a many-body calculation [9MTl1| . However, even if the 
band-gap had been corrected, the calculation of the opti- 
cal properties is not easy. The standard approach based 
on the adiabatic local density approximation ( ALDA) 0] , 
for example, dramatically underestimates the low energy 
peak - commonly referred to as the "excitonic peak" - 
in the optical spectrum. Improvements on the ALDA 
such as the adiabatic extension of the GGA do not fare 
much better. The exact exchange approach 01 has been 



found to catch the excitonic effect in silicon (12|, but 
this was achieved at the cost of artificially restricting the 
set of states included in the calculation to avoid the "col- 
lapse" of the spectra. So far, the most consistent ab-initio 
scheme leading to results in good agreement with exper- 
iment has been the recasting of Bethe-Salpeter equation 
as an equation for a two-point function within the frame- 
work of TDDFT but even this approach remains 
computationally very demanding. 

In recent years, a new class of approximate functionals 
has emerged in ground-state DFT. These are known as 
meta-GGA (MGGA) functionals and their defining char- 
acteristic is to depend not only on the density and its 
gradient, but also o n the non-interacting kinetic energy 
density r(r) 



13Hl6l|. At first sight, the dependence on 



r(r) seems to contradict the general statement that the 
xc potential is a functional of the density. But, it must 
be kept in mind that -r(r) is determined by Kohn-Sham 
orbitals which, in turn, are nonlocal functionals of the 
density. Thus the MGGA functionals are still function- 
als of the density, but intrinsically nonlocal ones. Their 
power stems entirely from this fact. 

In this Letter we show that an adiabatic approxima- 
tion based on meta-GGA functionals leads to very sig- 
nificant improvement in the calculation of optical prop- 
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erties. The fact that meta-GGA functionals can lead to 
improvements in the calculation of the KS band gap has 
been known for some time |16j |. What we add here to 
that knowledge is the realization that these functionals 
can also produce accurate optical spectra. And since 
the use of the adiabatic approximation automatically ex- 
cludes retardation effects, we conclude that the primary 
reason for the success of the meta-GGA functionals is 
the improved treatment of the long-rangedness in the xc 
potential. This long-rangedness (often referred to as "ul- 
tranonlocality" ) has long been known to be a problem 
in TDDFT, especially so in the applications to extended 
systems. While its strength could be inferred from fits to 
experimental spectra [H , none of the approximations de- 
veloped so far could deal with it satisfactorily. We believe 
that the use of meta-GGA functionals is a breakthrough 
in the handling of ultranonlocality and paves the way to 
efficient first-principle calculations of the optical proper- 
ties of semiconductors and more complex materials. 

Formulation - The crucial quantity targeted in 
TDDFT is the density-density response function 
x(r,r',u;), which is related to the non-interacting KS re- 
sponse function x s (r, v' ,u>) by the equation Q 



defined as the functional derivative of the dynamic xc 
potential v xc {r,u>) with respect to the dynamic particle- 
density. In order to calculate f xc we start from the ex- 
pression for the xc energy within MGGA as 



E xc = I e xc [n(r),Vn(r),r(r)]dr, (2) 



where the xc energy density e xc is a local function of its 
three arguments, n(r) is the particle density, 



r(r) = i]T/ Q |V^(r)| 



= ^/ Q e Q |V Q (r)| 2 - « a (r)n(r) + -V 2 n(r) (3) 

a 

is the non-interacting kinetic energy density, ipa, £a, and 
f a are the KS orbitals, their eigenenergies, and the occu- 
pation numbers, respectively, and v s (r) is the static KS 
potential. The second equality in Eq. follows from 
the KS equation. With the use of the definitions of the 
xc potential and the xc kernel as the first and the second 
functional derivatives of E xc with respect to density, we 
derive from Eq. ([2]) 



X 1 (r,r',u) = Xs 1 {r,r',u)-f xc (r,r',w)- 



(1) 



where f xc (r,r' ,u) = 6v xc (r,ui)/5n(r' ,u>) is the xc kernel, 



/ x dt xc , . „de xc , . f de xc , ,,<5r(r') , . ... 
v xc (r) = ^ (r) -V^(r) + y -^(r')^dr', (4) 
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(5) 



The xc potential of Eq. (j4]) has been thoroughly ad- 
dressed in Ref. Il7| and our focus will be the xc kernel 
of Eq. ([5]). With the use of the standard perturbation 
theory, the functional derivatives of r evaluate to fl8j 



5t(t) 
Sn{r>) 



-v s (r)5(r - r') 

J ff(r,r")x 8 - 1 (r",r')dr" + iv 2 <5(r-r'), (6) 



5 2 t(y) 



8n{v')5n{v") 



+2 Jf(v, n ,r 2 ) X ; 1 (n Z)*; 1 (r 2 ,r")dndr2 
+ J ff(r,r 1 )Xrt 1 (ri,r',r ,, )dri J (7) 



where 



2 ^ e «- £ /3 



(8) 



F(r,ri,r 2 ) 



[C(r 2 )^ 7 (r 2 ) 

x ^(r)^(r)V ) a(i"i)V'K r i) + ( r ^ r i)] 
- E ull\& [l^(n)| 2 C(r 2 )^(r 2 ) 



x Va(r)^(r) +(ron) + (n o r 2 )] , (9) 

and \~2 in Eq. (0 is the inverse of the quadratic KS 
density-response function x S 2(r,r',r") = Sv J r "sv t ( T ") - 
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Equations ©-(HI) together with the explicit KS re- 
sponse functions [l8[ constitute the complete solution to 
the MGGA-based xc kernel f xc (in the adiabatic approx- 
imation) . 

Ultranonlocality - In reciprocal space, the xc kernel 
becomes a matrix in the reciprocal vectors fxc,GG' (l) • 
Whether or not the MGGA for f xc provides an improve- 
ment over conventional approximations depends on the 
presence or absence, in the optical limit q — > 0, of a 
singularity of the type / 2 - c ,oo(q) — a/q 2 [sj]. Obviously, 
LDA and GGA [the first three terms in Eq. §5§] do not 
have such a singularity. With the neglect of the non- 
singular terms in the matrix form of Eq. ([5]) and without 
the "local-field effects" , we simplify /,„,, to [H| 



icGG 



(q) 



de xc 
' dr 



X.GG'(q). 



(10) 



where the overlinc denotes the average over the unit cell. 
The right-hand side of Eq. (fTU I contains the singularity 
in question because xJ 1 does 19j. Focusing on the 00 
component, we finally get 



8e 



(11) 



Considering that de xc /dr is almost the same for Si and 
Ge, neglecting for a moment the local-field effects, and 
neglecting the unity compared to the static dielectric 
function of a semiconductor, we see that Eq. (fTTj) is in 
agreement with the empirical rule of a being inverse- 
proportional to the dielectric function [ll| . We also note 



that the ultranonlocality we find seems to be the first ex- 
plicit demonstration of the fact that the kinetic energy- 
dependent functionals are not in practice semi-local in 
the density Q. 

Choice of functionals and calculation of optical proper- 
ties - Having established on the fundamental level that 
the adiabatic meta-GGA-based TDDFT does account for 
the ultra-nonlocality in crystals, we now turn to numer- 
ical calculations. First we note that onl y t he group of 
functionals that provide E xc (e.g, VS98 [H and TPSS 
14j) rather than those providing v xc directly (e.g., BJ06 
15| and TB09 (HI) can be used to build f xc , since for 
the functionals of the latter group the corresponding E xc 
does not exist [ltij . We have used two well established 
MGGA functionals VS98 fij and TPSS for the cal- 



culation of both the ground-state with v xc of Eq. (UJ and 
f xc of Eq. ([TU]) [20| . The resulting values of the key quan- 

are listed in Table |U At 



tity de xc /dr entering Eq. ([TO]) 
first glance surprisingly, the values found with the two 
different functionals differ drastically: The r-dependent 
part of the TPSS functional was found negligible every- 
where over the unit cell. In the Supplementary material 
[lij . we analyze the r-dependence of VS98 and TPSS 
functionals to the conclusion that for the latter it is very 
weak. Accordingly, we argue that while well tuned to 



TABLE I. The average over the unit cell of the derivative of 
the exchange, correlation, and xc energy density with respect 
to the kinetic energy density found for Si and Ge with the 
VS98 GJl and TPSS El functionals. 



VS98 


TPSS 


de x /dr 


de c /dr de xc /dr 


de x /dr 


de c /dr 


de xc /dr 


Si 0.122 


-0.226 -0.104 


4.60xl0 -3 


-1.15X10 -4 


4.49X10 -3 


Ge 0.135 


-0.241 -0.106 


2.94xl0 -3 


8.73xl0" 5 


3.03X10" 3 



yield accurate E xc , TPSS performs unsatisfactorily with 
respect to its r-derivative. A clear reason for the weak 
r-dependence of TPSS can then be easily identified: This 
functional is tuned to (i) the nearly free electron gas 
(NFEG) and (ii) the one and two electron systems |14j . 
In both cases, due to the gradient expansion of the kinetic 
energy of NFEG and to the von Wcizsackcr's formula for 
the kinetic energy of one and two electron systems, re- 
spectively, t is (semi)-local in density, which leads to the 
local theory with respect to f xc and zero a (Cf. [H|). 
On the other hand, the VS98 functional is designed to 
work better in the strong rather than the weak inhomo- 
geneity case [13(, which qualitatively explains its success 
in yielding realistic values of a. Accordingly, we use the 
latter MGGA functional in our calculations below. 

We calculated the KS band-structure and the micro- 
scopic density-response matrix of Si and Ge with the 
full-potential linear augmented plane-wave (FP-LAPW) 
method and the VS98 MGGA xc functional [Hj]. The 
supporting results for zincblende semiconductors are pre- 
sented in |18| . The real and imaginary parts of the macro- 
scopic (q = 0) dielectric function are presented in Figs.Q] 
and [21 It is evident that the inclusion of the the non- 
local f xc of Eq. (fTOj) via the MGGA greatly improves 
the agreement between the theory and experiment, in 
particular, making the excitonic peak considerably more 
pronounced. 

It is instructive to draw a parallel between our ap- 
proach and that of Ref. 0. Before the inclusion of f xc , 
both methods produce single-particle spectra that un- 
derestimate the intensity of the excitonic peak. Then, 
with the inclusion of the many-body interactions through 
f xc , the spectra are red-shifted and the excitonic feature 
grows. The fundamental differences between the two ap- 
proaches are that: (i) We remain all the time within the 
framework of TDDFT, while Ref. © uses a combination 
of TDDFT with the GW approximation of many-body 
theory; (ii) While the quantity a in Ref.© was a fitting 
parameter, we have for it an explicit expression, Eq. (|lip . 
Moreover, for Si, Eq. ([IT]) evaluates to a = —0.267, which 
compares reasonably well with the best fit value of Ref. © 
of a= -0.2 (24 1 . 

In conclusion, we have developed the adiabatic 
TDDFT formalism for the kinetic energy dependent 
(MGGA) exchange-correlation functionals. In con- 
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FIG. 1. (color online) Dielectric function of silicon. Thin solid 
(red online) line is the result obtained with MGGA band- 
structure and including the many-body interactions through 
f xc of Eq. pop . Dashed (green online) line is the result ob- 
tained with MGGA band-structure but with f xc = (RPA). 
Dotted (blue online) line is obtained with LDA band-structure 
within RPA. Thick solid line is the experiment from Ref. I2H , 
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FIG. 2. (color online) The same as Fig. [1] but for germanium. 



trast to LDA and GGA approximations, the result- 
ing exchange-correlation kernel f xc is shown to exhibit 
the singularity of the type a/q 2 , which is a neces- 
sary feature for a theory to describe the excitonic ef- 
fect in crystals. Our calculations performed for a num- 
ber of the diamond-structure and zincblcndc semiconduc- 
tors demonstrate the high promise of the MGGA-bascd 
exchange-correlation functionals as a new tool in the ar- 
senal of TDDFT methods. 
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AUXILIARY MATERIAL 



TO THE PAPER BY V. U. NAZAROV AND G. VIGNALE 

"OPTICS OF SEMICONDUCTORS FROM THE META-GGA-BASED TIME-DEPENDENT 

DENSITY-FUNCTIONAL THEORY" 

1. Derivation of Eqs. ©-([9]) 

To evaluate the functional derivatives of r explicitly, we note that the standard perturbation theory yields 

Se n , , , M2 



Sv s (r) 



IV> Q (r)| 2 , (A.l) 



M r f^ a z a -ep 



and by virtue of the functional chain rule 

dn(r')~J Sv s (r") 5n(V) ~J 5v s (r") Xs 1 ' ' ' 1 ' 

from Eq. using Eqs. (|A.1[) and (|A.2|) . we straightforwardly (though, with a lengthy algebra) arrive at Eqs. © 
and (0. Wc also mention the properties 

J H(r,n)dr = J ff(r,n)<2ri = 0, (A.4) 
y > F(r,ri,r 2 )dr = x.7 1 (i'i,r2), (A.5) 

which directly follow from the definitions ([5]) and © with account of the orthogonality of the KS orbitals, and which 
will be used below. 

2. Explicit forms of the density-response functions 
In equation (|7|), x^2 ^ s ^ ne inverse of the quadratic KS density-response function. A convenient relation holds 

X~l(r,r',r") = - J xj^r, r 3 )x s2 (r 3 , r 4 , r 5 )x7 1 (r 4 , r')x7 1 (r 5 , r")dr 3 dr i dr 5} (A.6) 

which can be easily proven by functional differentiation of the identity XsXj 1 = 1- Explicitly, in terms of the KS 
orbitals, eigenenergies, and the occupation numbers 

Xs (r, r') = £ ^-Ilr Q (r)Mr)Mr'W P (r r ), (A.7) 

Xs2 (r,r',r") = £ /"J 7 * [|^(r")| 2 " l<MO| 2 ] CW^(r)^(r')^(r') + (r O r") + (r' O r") 

+ E 7 IelzIJL ^(r")^ 7 (r")^(r)^(r)^ a (r')^(r') + (r r") + (r o r') + (r -> r" -> r' r). (A.8) 
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3. Derivation of Eq. (flO]> 

Although Eq. (|10p can be obtained directly as an approximation to Eq. (|5]), it is more transparent to take a step 
back to Eq. (j4]). Neglecting all the Fourier coefficients of the function (r) but the zeroth (which is supported by 
our calculations for Si and Ge) we have from Eq. (0} 



de. 



v xc r « — — r - V— — (r, , . 

an avn or y on(r) 



Using Eq. ([6]) and the property (|A.4|) of the function H of Eq. flSJ, we can write 



<Jn(r 



-dr' 



-u s (r). 



Another (instructive) way to prove Eq. (jA.101) is to note that by the minimum principle 



_6E_ 
6n(r) 







8n(r) 



J T (r')dr' 



+ v ext (r) + v H (r) + v xc (r), 



(A.9) 



(A.10) 



(A.11) 



where E is the energy of the system. Equation (|A.llj) immediately yields Eq. (|A.10|) . 

Substituting Eq. (|A.10[) into Eq. (|A.9[) and taking the functional derivative of v xc (r) with respect to the density, 
we arrive at 



Sn(r') I dn 



(r) - v aw (r) 



-v s (r) 



8 de xc de xc 8v s (r) 



Sn(r') dr dr Sn(r') 



(A.12) 



Finally, with the use of Eq. (J6j), it is straightforward to show that all the terms in Eq. (|A.12|) but the last are not 
singular as 1 /q 2 in the reciprocal space. By neglecting those terms and recalling the relation 



-1/ / \ Sv a (r) 



(A.13) 



we conclude the derivation of Eq. (|T0 



4. Comparison between VS98 and TPSS 



Figure fA . 1 1 demonstrates that while producing almost the same xc energy density in crystalline silicon (lower panel), 
VS98 and TPSS functionals do it with very different 'weights' of the dependence on n and Vn on one hand, and on 
r, on the other. In particular, TPSS functional exhibits almost no dependence on r in the range of parameters of 
this system (upper panel), hence, doing its job of approximating e xc by means of n and Vn only. In other words, 
TPSS functional is almost a GGA one within this range of parameters. On the contrary, VS98 functional depends on 
t crucially, hence, being really nonlocal with respect to n. 
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5. Results for more semiconductors 



Using Eq. (|T0|) . we have also performed calculations for zincblende semiconductors with the results shown in Fig. IA.2l 
To facilitate the calculations and to include f xc on top of the as accurate as possible fundamental gaps, in the ground- 
state calculations we have used the highly efficient and accurate Tran and Blaha xc potential (TB09) while f xc 
has been included with VS98 as before. In order to demonstrate that the differences with the consistently used VS98 
throughout (Figs. [T] and [2]) are minimal, we also present results for Si and Ge obtained by this slightly simplified 
method. We note that TB09 MGGA xc potential cannot be used for constructing f xc since it gives the xc potential 
directly which docs not correspond to any xc energy-density fTpj ] . 
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FIG. A. 2. Dielectric function of Si, Ge, GaAs, GaP, and InAs semiconductors. Red line is the result obtained with the 
inclusion of the many-body interactions through f xc of Eq. (|10|) using the VS98 MGGA xc functional while the ground-state 
band structure was obtained with TB09 xc potential. Blue line is obtained with LDA band-structure within RPA. Thick solid 
line is the experiment from Ref. |23| . 



